Equation of State, Stability, Anisotropy and Nonlinear Elasticity of Diamond-Cubic 
(ZB) Silicon by Phonon Imaging at High Pressure 
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Experimental phonon imaging in diamond anvils cell is demonstrated to be an adequate tool to 
extract the complete set of elastic constants of single- crystalline silicon up to the ZB— >■ /?— Sn transi- 
tion (10 GPa). Contrary to what was commonly admitted, we demonstrate that the development of 
the strain-energy density in terms of strains cannot be stopped, for silicon, after the terms contain- 
ing the third-order elastic constants. Nonlinear elasticity, degree of anisotropy and pressure-induced 
mechanical stability of the cubic silicon structure are thus revisited and investigated in more detail. 

PACS numbers: 62.20.de, 62.50.-p, 64.30.Jk 



I. INTRODUCTION 

The determination of interatomic forces in solids with 
respect to atomic displacements is one of the most gentle 
and direct way to probe the short-range repulsive poten- 
tial. Experimental study of elastic properties under high 
pressure has thus become an essential part of solid-state- 
physics, although such measurements are still a nontrivial 
task. For example, the measurement of elastic constants 
as a function of the pressure (i.e. the crystal volume) can 
serve to derived a precise and hypothesis- free equation of 
state (p,V,T). It is also useful to understand structural 
and electronic instabilities in solids. Last but not least, 
the well-known accuracy of ultrasonics measurements has 
been shown to allow the determination of both second- 
and higher-order elastic constants, giving rise to the char- 
acterization of nonlinear and anharmonic properties of 
solids. 

Already studied for more than fifty years through the 
pressure variation of sound waves under moderated pres- 
sures (less than 1 GPa) , the determination of nonlinear 
properties has recently attracted great interest, in con- 
jonction with the development of nanoelectromechani- 
cal and nanoelectronic devices. In such systems, the 
large surface contribution to the total energy modifies 
the properties and, in particular, raises the intrinsic non- 
linear properties. In a recent paper, Lopuszyirski et a^ 
have shown through DFT calculations that the knowl- 
edge of pressure derivatives of elastic constants is a clear 
prerequisite to a complete explanation of semiconductors 
nanostructures properties. Tang et al^ came to the same 
conclusion through a calculation of lattice dynamics in 
silicon nanostructure, one of the most interesting system 
for applications as sensors or communication technolo- 
gies. 

Extrapolating the pioneering work of McSkimin and 
Andreatch in 1964 on the elastic properties of cubic 
silicon^ and germanium^ up to 0.2 GPa, the hydro- 
static pressure dependence of the 3 independent Second- 
Order Elastic Constants (SOEC) of covalent semicon- 



ductors has always been considered to be linear. Rig- 
orously, these experimental results only give an access to 
the determination of few Third-Order Elastic Constants 
CijK (TOEC) combinations, constants which can sim- 
ply be expressible in terms of SOEC pressure derivatives 
dC/j/dp. To determine the six independent TOEC of a 
cubic crystal, one also required experiments under uni- 
axial stress. However, difficulties inherent to the non- 
hydrostatic technique^ are known to produce doubtful 
conclusions^ and/or severe disagreements between au- 
thors. Even more conflicting are the inconsistencies be- 
tween hydrostatic and uniaxial data in a same work^i. To 
avoid any controversies, we here restrict our attention to 
the nonlinear effects of a single crystal of silicon under 
hydrostatic pressures. 

Reaching higher pressure as possible is crucial since the 
rather subtle contribution of high-order elastic constants 
is known to increase with increasing pressure. From an 
experimental point of view, our main goal was thus to 
improve our knowledge on the elasticity of cubic silicon 
through an highly accurate sonar measurement over the 
largest pressure domain of observation. To follow this 
objective, we have used the purpose-built technique of 
picosecond acoustics phonon imaging at high pressure. 
This method, described in the first part of the article, 
combines the great accuracy of ultrasonic laser technique 
with the capabilities of diamond anvils cell (DAG). It al- 
lows the measurement of the complete set of elastic con- 
stants up to the transition pressure (10 GPa) of a defect- 
free and undoped single-crystalline silicon sample^. In a 
second part, we show that such measurements may serve 
as a starting point to discuss the validity of the linear re- 
gression hypothesis of C/j(p), commonly admitted. Non- 
linearity and anisotropy of the Si elastic properties as a 
function of pressure are also discussed in terms of phase 
stability. This may be used as a stringent test of accuracy 
(or even validity) for state-of-the-art simulations, or even 
as a significant information to improve the modeling of 
silicon-based nanoelectromechanical and nanoelectronic 
devices. 



II. PHONON IMAGING IN DAC : 
ARGUMENTS AND PRINCIPLE 

The main reason why acoustic measurements on sin- 
gle crystal at high pressure are quite sparse in litera- 
ture is technical. Traditionally, the frequency of sound 
waves generated by ultrasonic transducers lies typically 
in the MHz range, so that any sonar-based measurements 
would require millimetric sample sizes. Such high di- 
mensions render the piezoelectric transducer technique to 
only be performed when combined with a large-volume 
celli^. This, consequently, limits the pressure domain. 
To overcome this limitation, an attempt to implement 
acoustic GHz interferometry^-'^ with DAC has been done. 
This promising method is however still limited in pres- 
sure due to the presence of diffraction effect or acoustic 
phase shifts at the diamond-sample interface under non- 
hydrostatic stressesiS.. The third experimental scheme is 
the well-known Brillouin scattering, an optical method 
that can easily be adapted to DAC^"^, but which is un- 
suitable in the case of silicon : in Brillouin scattering 
method, transparent sample can exclusively be studied. 

More recently, picosecond acousticsi^"— has been de- 
veloped to enable highly accurate mechanical properties 
of opaque materials as metal liquidsi^, polycrystals;^^, or 
single-crystali2 in DAC. In this method, the measure- 
ment is similar to the pulse-echo ultrasonic technique 
(travel time determination) with the advantage of the 
optical methods (no contact ; no bonding effect). The 
principle of picosecond acoustics in DAC, described in 
referencesiiiS, is based on the absorption of the light 
pump pulse which sets up a local thermal stress near the 
surface. In the acoustic far field diffraction limit (i.e. the 
laser illuminates the sample with a spot size much higher 
than the film thickness) , only longitudinal acoustic strain 
field is created. In order to circumvent the "shear less" 
problem, we further focus our attention on the poten- 
tiality of picosecond ultrasonics to generate and detect 
shear waves in sample embedded in a DAC. Keynote is 
to generate lateral compressive stresses (producing inter- 
nal diffraction) through a minimization of the source area 
with respect to the characteristic acoustic wavelength. 

Original development of acoustics using small source 
and point detector have revealed some striking proper- 
ties of phonons in anisotropic materials, such as pho non- 
focusing effect—. In crystals, the direction of acoustic 
energy flow does generally not coincide with the wave 
vector direction k. Consequently, a source producing 
an uniform angular distribution of k will generate an 
anisotropic propagation of elastic energy. Considering 
a crystal with plane parallel surfaces, a point acoustic 
source at one surface will produce non-uniform focus- 
ing pattern at the opposite face. This produces pat- 
terns related to the complex topology of the group ve- 
locity surface (also called wave surface). Such intrigu- 
ing imaging can however be calculated from elastic the- 
ory (within the continuum model) taking into account 
the relation between directions of k (or phase velocity 



V — a;(k)/k) and directions of energy flow (or group ve- 
locity v*^ = 9a;(k)/9k). Experimentally, the earlier im- 
ages of flux energy pattern, obtained in heat-pulse exper- 
imenta^i^ have shown that a massive amount of data may 
be extracted from successively recording acoustic wave- 
front snapshots at different times. In recent years, the 
technique of picosecond laser acoustics has been shown 
to be well suited for such purpose : femtosecond dura- 
tion of the laser pulse and the possibility to focalize the 
beams onto the surface plane of the sample give rise to an 
outstanding resolution in space and time. In particular, 
the measurements of group velocities using picosecond 
acoustics has been shown to allow the determination of 
the complete set of elastic constants of anisotropic sample 
at ambient pressure and high temperature^-. Following 
this approach, and taking into account the considerable 
progress made the last five years in ultrafast acoustics^^, 
we have here implemented the wavefront imaging method 
in the case of picosecond acoustics in DAC. 



III. EXPERIMENTAL SET-UP 

Ultrashort pulses of 100 fs are generated every 12.6 ns 
by a Ti:Sapphire laser (800 nm). The laser beam is split 
into pump and probe beams. The pump is focused on 
one surface of the sample, whereas the probe is focused 
on the opposite one. As soon as the pump laser pulse 
reaches the surface, it creates a sudden and small tem- 
perature rise (of about 1 K). The corresponding thermal 
stress generated by thermal expansion relaxes by launch- 
ing acoustic strain fields. After propagation along the 
sample, both thermal and acoustic effects alter the op- 
tical reflectivity of the sample in two ways: the photo- 
elastic effect, and the surface displacement (as the acous- 
tic echo reaches the surface). The first modification con- 
tributes to the change of both real and imaginary parts of 
the reflectance, whereas the second one only modifies the 
imaginary contribution. In a pure thermo-elastic model, 
the time and space refiectivity change Ar{t) can be repre- 
sented as a function of the photo-elastic coefficient dn/dz 
and the surface displacement Mo(i) as : 
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where n is the optical index of the sample. 
The variation of refiectivity as a function of time is de- 
tected through the intensity modification of the probe 
delayed from the pump with a different optical path 
length. The uncertainty on the absolute position of the 
delay line being less than 1 /xm, the determination of 
the pump/probe time delay is better than 20 fs. The 
detection is carried out by a stabilized Michelson inter- 
ferometer which allows the determination of the refiectiv- 
ity imaginary part change^^. Microscope objectives are 
here necessary in order to focus both pump and probe 



beams down to 3 /zm (objectives with a typical working 
distance of about 20 mm and an optical aperture of 0.42, 
well adapted to the DAC environment, have been used). 
Finally, an image is obtained by scanning the reflectiv- 
ity of the whole sample surface through the displacement 
of the probe objective in the plane perpendicular to the 
beam using a X-Y piezoelectric stage (giving rise to im- 
ages of 100 /im X 100 /im, with an accuracy on the ab- 
solute position better than 1 ^m). We emphasize that, 
albeit the experimental set-up is similar as the one used 
in a sonar configuration (see Ref^^-), the surface displace- 
ment of the sample is here determined by scanning the 
surface at fixed pump-probe delay (i.e. at a given time 
along the acoustic propagation). 

Extracted from a large single-crystal of silicon, a thin 
platelet oriented along [100] with surface of about 70*70 
/im^ was mechanically polished. A thin film (50 nm) of Al 
was sputtered on both sides to serve as transducers. This 
sample was loaded into the experimental volume of the 
DAC. Neon was used as pressure transmitting medium. 
Using the well known longitudinal velocity of Si along 
the [100] direction at ambient conditions^, the thickness 
of the crystal platelet was measured to be 42.2(1) ^m 
from picosecond measurement in the classical configura- 
tion (i. e. scanning the delay at fixed pump-probe rel- 
ative position). Pressure was classically measured using 
the fluorescence emission of a 5 /xm ruby sphere^'^ placed 
close to the sample in the gasket hole. The accuracy was 
better than 0.1 GPa at the maximum pressure reached. 



IV. RESULTS AND DISCUSSION 
A. Elasticity at ambient conditions 

Using the density p=2.331 g.cm^"^ and a set of initial 
elastic constants values, the simulated slices of the wave 
surface in the (100) plane for silicon at ambient condi- 
tions are obtained in two steps. Firstly, the Christoffel 
equation is solved for a set of wave vectors k lying within 
45 degrees (cubic symmetry) around the [100] crystallo- 
graphic direction. The slowness curves are then gener- 
ated for the three acoustical polarizations and, in a sec- 
ond step, used to calculate from ray theory2£ the wave 
front curves within a surface cut in the plane (100). In 
Figure [H a typical 3-D phonon imaging pattern for the 
three acoustic polarizations is given. 

The experimental and calculated intercepts between 
3-D wave front and (100) plane of silicon at ambient con- 
ditions is shown in Figure [H The simulated process of 
the wave surfaces, where the elastic constants are the 
fitting parameters, well reproduces the experimental pat- 
tern and gives rise to an elastic tensor in very good agree- 
ment with previously published data^. We here empha- 
size that, for a given thermodynamical condition (here 
ambient), different patterns observed for different pump- 
probe delays can be used to renew the fitting process in 
order to determine the complete set of elastic constants 




FIG. 1: Calculated group velocity surfaces near the [100] 
symetry axis of ZB-cubic silicon. Red, blue and green dashed 
lines correspond respectively to the longitudinal, fast and slow 
transversal group velocities. 



with a higher accuracy. 



hjflkkL 







FIG. 2: Experimental phonon imaging pattern in the (100) 
plane of silicon at ambient conditions (pump-probe delay time 
of 0.4 ns). Red, blue and green dashed lines correspond re- 
spectively to the longitudinal, fast and slow transversal group 
velocities using Cii= 165.7 GPa, Ci2= 63.9 GPa and C44= 
79.5 GPa. 

Note that the fast transversal mode was not experi- 
mentally detectable and thus not been taken into account 
in the elastic constants fitting process. The absence of 
this mode (polarized into the (100) plane) can be eas- 
ily understood using Eq. [TJ where the imaginary part is 
mainly dominated by the surface displacement perpen- 
dicular to the surface. Consequently, using our interfero- 
metric configuration, any surface displacement could be 
detected except when it is perpendicular to the surface. 
A calculation of a [100] projection of polarization vec- 
tors, shown in Figure [31 illustrates the absence of the 



fast transversal mode for which the corresponding dis- 
placement is almost zero. 




FIG. 3: Calculation of polarization vectors (eigenvectors of 
the Christoffel equation) projected along the [100] direction, 
perpendicular to the surface where the probe is focalized. 
Red, blue and green dashed lines correspond respectively to 
the longitudinal, fast and slow transversal polarization. 



B. Nonlinear ity and anisotropy at high pressure 

At high pressure, a step forward of the previous sec- 
tion, typical experimental and calculated slices of the 
wave surface in the (100) plane are given in Figure 01 




FIG. 4: Top : experimental phonon imaging patterns in the 
(100) plane of silicon at 7.75 GPa at two different probe-pump 
delays. Bottom : same as top with superimposed calculation 
curves for longitudinal, fast and slow transversal group ve- 
locities (red, blue and green dashed lines respectively) using 
Cii= 196.9 GPa, Ci2= 104 GPa and C44= 80 GPa. 



Taking into account the pressure dependence of the 
density, we have here used a classical procedure to deter- 
mine the thickness of the crystal plate at each pressure. 
The basic idea can be briefly described as follows. As- 
suming that length and density are known at a given 
pressure p (ambient pressure, for example), all velocities 
and elastic moduli can be directly deduced from the pi- 
cosecond measurements. Using p- values of density and 
thickness, and p-|-Ap wavefront snapshots for different 
pump-probe delay, the length and the density of the sam- 
ple at p-|-Ap can be computed. These values are then 
used to deduce a first approximation of the elastic mod- 
uli at the new pressure p-|-Ap. The values of the latter 
are finally used as starting points in an iterative process, 
until convergence is reached. This process is quite robust 
since the variation of the sound velocity is mainly due to 
that of the elastic moduli (i.e. its dependence on length 
and density is only of second order) . 

Figures [S] show elastic constants of cubic diamond (ZB) 
silicon as a function of pressure up to the ZB-^- /3— Sn 
transition (10 GPa). An excellent agreement is observed 
with computations results using DFT-LDA total-energy 
method^-. The comparison is also very good with tight- 
binding simulations^^, except for C44 which calculated 
ambient pressure value is far from the well established 
value of 80 GPaSS. Whereas this discrepancy may be 
due to an inadequate model parametrization, a good 
agreement between theoretical and experimental pres- 
sure derivative of C44 is nevertheless observed. Present 
data and 1960's ultrasonics results^ are also in excellent 
agreement (in the pressure domain where the compari- 
son can be done, say less than 0.2 GPa). In this pres- 
sure range, the infinitesimal strain theory does not hold 
anymore and, within the accuracy of the measurements, 
second-order moduli C/j are observed to vary linearly 
with pressure. However, a comparison between our data 
at higher pressure (typically more than 3 GPa) and a 
linear extrapolation of the McSkimmin^ results is worse. 
Beyond the classical quasiharmonic approximation, the 
pressure dependence of the third-order moduli Cijk of 
silicon is here needed to interpret such disagreement. A 
second-order polynomial regression of present data yields 
to (data in GPa) : 

Cii(p)=165.7-l-4.73p-0.09p2; 

Ci2(p)=63.60+5.78p-0.12p2; 

C44 is however unaffected by a variation of pressure, ly- 
ing around 80 GPa within the experimental uncertainty. 

Consequently, the pressure change of three inde- 
pendent relations for Cijk (here labeled Ei{p)) is 
calculated— to be (data in GPa) : 

Ei(p)=Ciii+2Cii2=-1860+28p; 

E2(p)=2Cii2+Ci23=-1120+31p; 

E3(p) = Ci44 + 2Ci66=-747. 

A good agreement is observed (see Table H]) between 
the zero-pressure values Ei(0) with previously published 
data. 

Whereas an evaluation of any of the six independent 
Cijk specifically is beyond the capability of elastic mea- 




TABLE I: Comparison of experimental (exp) and calculated 
(calc) relations Ei for Cuk (in GPa) of silicon at ambient 
conditions. 
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FIG. 5: Pressure dependence of cubic silicon elastic constants 
Gil (up), G44 (middle) and G12 (down) at 300 K. The data for 
Gil corresponds to experiment of phonon imaging (open cir- 
cles) completed by an experiment in a classical high pressure 
picosecond set-up^. Present results (with uncertainty given 
by the symbol size) are compared with previous data (blue 
solid line : Tight-binding DFT calculations^^, green solid line 
: LDA-DFT calculations^^, and dot black line : extrapolation 
of low-pressure ultrasonics data*). The red solid line corre- 
sponds to a second-order polynomial fitting of the phonon 
imaging data. 
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clearly point out the contribution of the fourth-order 
elastic constants (FOEC) Cijkl- In qualitative agree- 
ment with our present observation, Prasad et aP^ have 
evaluated Cnn to he around 12 TPa, high value which 
entails the FOEC contribution to be substantial at high 
pressure. This clearly demonstrates the need to account 
for high-order elastic constants in the study of the anhar- 
nionic properties of silicon, of most importance for sim- 
ulations (test of validity for explicit crystal potential), 
or for applied physics (harmonics generation and distor- 
tion effects) as sound waves propagating in silicon-based 
devices. 

As well as nonlinearity, elastic-anisotropy is also known 
to shed light on the mechanical properties of solids as dis- 
location dynamics, plastic deformations or structural sta- 
bility. For cubic-symmetry as silicon, the 1930's Zener— 
definition of the acoustical anisotropy A=2Cii/{Qii- 
C12) has been recently revisited^^ in order to quantify the 
single crystal elastic anisotropy from an universal point of 
view A^ = 6/5(v]4— 1/vA)^ (note that for an isotropic 
crystal, A^l but A^=0). As illustrated in Fig. H the 
degree of elastic anisotropy of silicon is increasing with 
increasing pressure. This result well agrees with ab-initio 
calculationa^li^ and could be related to the observation 
that the transverse phonon frequency (i.e. the gruneisen 
parameter) has been found to be the unique mode that 
decreases with pressure^^. 



C. Equation of state and structural stability 

For cubic crystals, the bulk modulus B is a simple func- 
tion of the elastic constants : 

B=(Cii+2Ci2)/3 

Using the experimental values of C/,/(p), highly accu- 
rate and free- hypothesis equation of state B(p) has been 
determined (see Fig. [T]). 

A second-order polynomial curve well reproduces the 
experimental pressure dependence of the bulk modulus 
with : 

B(p)=97.60+5.74p-0.16p2 

This result in excellent agreement with a fourth-order 
regression of the elastic energy with respect to the strain 
giving : 

B'=dB/dp=-(Ei+2E2)/(9B)=5.08-0.24p 

Lattice stability and sequence of phase transformations 
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FIG. 6: Universal elastic anisotropy factor A of ZB silicon as 
a function of the pressure (open circles) compared with previ- 
ous data (blue solid line : Tight-binding DFT calculations"'*, 
green solid line : LDA-DFT calculations^i, and dot black line 
: extrapolation of low-pressure ultrasonics data^). 



FIG. 8: Elastic stability criteria B2 as a function of pres- 
sure. Extrapolation of present results (black line) are com- 
pared with previously published calculations (blue solid line : 
molecular dynamics simulations^^ and red solid line : LDA- 
DFT calculations^! . 
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FIG. 7: Pressure dependence of bulk moduli in single- 
crystalline ZB silicon. The uncertainties in pressure and 
bulk modulus are within the symbol size. Present data are 
compared with previously published results (blue solid line 
: Tight-binding DFT calculations-^, green solid line : LDA- 
DFT calculations^^ , and dot black line : extrapolation of low- 
pressure ultrasonics datai). 



in semiconductors can be understood in the framework 
of phonons instabilities'^^, which for purely covalent com- 
pounds as silicon can be expressed in terms of mechanical 
stability of the stressed lattice. In a cubic crystal, the 
condition of a positive density of elastic energy trans- 
forms to the following criteria : 

Bi=Cii+2Ci2 >0; 

B2=Cii-Ci2 >0; 

B3 = C44 >0. 

Bi is directly connected to the bulk modulus and does 
not allow to discuss the phase stability in terms of sym- 
metry change from cubic (ZB) to tetragonal (/3— Sn). 



However, B2 and B3 are referred to bulk shear and tetrag- 
onal shear moduli respectively. 

Theoretical studies of ZB structural stability of silicon 
have been undertaken through DFT— and molecular dy- 
namics simulations^^. In both works, the determination 
of the vibrational distortions at high pressure gives a crit- 
ical pressure (quite higher than the physical transition 
one) at which the ideal lattice become unstable against 
homogeneous tetragonal shear deformation, say B2. In 
our case, a tentative of extrapolation of experimental 
results can only be done for B2 since the pressure de- 
pendence of C44 was not enough accurate to be decently 
interpolate by other than a constant value. As can be 
seen in Fig. |S1 the violation of B2 occurs at around 120 
GPa, in good agrement with theoretical calculations. 

This result can be interpreted as the following : silicon, 
as mainly other tetrahedrally bonded semiconductors^^, 
is mainly stabilized by noncentral covalent interactions. 
From that point of view, the ZB— >■ /?— Sn path transi- 
tion has intrinsically a martensitic nature, ideal trans- 
formation that can not be effectively observed due to 
kinetic contributions. Consequently, the deformation of 
ZB-silicon under hydrostatic pressure can been seen to be 
driven by a tetragonal shear strain, but without reaching 
the structural instability. 



V. CONCLUSION 

We have presented a detailed experimental study of 
elastic constants of cubic ZB silicon as a function of hy- 
drostatic pressure. Phonon imaging in diamond anvils 
cell, an original experimental purpose-built set-up, al- 
lows to determine the subtle effect of high-order elastic 
constants and demonstrates the failure of the linear re- 
gression approximation dC/,7/dp=constant. From this 



study, nonlinear acoustics, anisotropy, equation of state 
as well as structural stability of a pure single-crystal of 
silicon have been extracted with an outstanding accuracy. 
So far, this technique is likely to have an impact on the 



study of dynamics of all single-crystals at high density. 
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